https://femiguez.github.io/nlraa-docs/nlraa-Oddi-LFMC.html
https://onlinelibrary.wiley.com/doi/pdfdirect/10.1002/ece3.5543
# install.packages("emmeans")
library(nlme)
library(nlraa)
library(ggplot2)
library(emmeans) # emmeans::contrast(), emmeans::emmeans()
library(tidyverse)
data(lfmc)
glimpse(lfmc)
Rows: 247
Columns: 6
$ leaf.type <fct> M. spinosum, M. spinosum, M. spinosum, S. bracteolactus, S. bracteolactus, S. bracteolactus, Gras…
$ time <int> 1, 1, 1, 1, 1, 1, 1, 1, 13, 13, 13, 13, 13, 13, 13, 13, 13, 26, 26, 26, 26, 26, 26, 26, 26, 26, 3…
$ plot <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ site <fct> P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, P…
$ lfmc <dbl> 238.4, 269.8, 325.2, 326.7, 257.3, 279.3, 91.9, 95.0, 215.8, 216.7, 205.3, 208.5, 225.7, 246.8, 8…
$ group <fct> SM plot 1, SM plot 1, SM plot 1, SS plot 1, SS plot 1, SS plot 1, GE plot 1, GE plot 1, SM plot 1…
lfmc
lfmc %>%
count(leaf.type)
lfmc %>%
count(time)
lfmc %>%
count(plot)
lfmc %>%
count(site)
lfmc %>%
count(group)
lfmc %>%
count(site, plot)
lfmc %>%
ggplot(aes(time, lfmc)) +
geom_point() +
facet_wrap(~group)
lfmc.gd <- groupedData(lfmc ~ time | group, data = lfmc, order.groups = FALSE)
lfmc.gd
Grouped Data: lfmc ~ time | group
leaf.type time plot site lfmc group
1 M. spinosum 1 1 P 238.4 SM plot 1
2 M. spinosum 1 1 P 269.8 SM plot 1
3 M. spinosum 1 1 P 325.2 SM plot 1
4 S. bracteolactus 1 1 P 326.7 SS plot 1
5 S. bracteolactus 1 1 P 257.3 SS plot 1
6 S. bracteolactus 1 1 P 279.3 SS plot 1
7 Grass E 1 1 P 91.9 GE plot 1
8 Grass E 1 1 P 95.0 GE plot 1
32 M. spinosum 13 1 P 215.8 SM plot 1
33 M. spinosum 13 1 P 216.7 SM plot 1
34 M. spinosum 13 1 P 205.3 SM plot 1
35 S. bracteolactus 13 1 P 208.5 SS plot 1
36 S. bracteolactus 13 1 P 225.7 SS plot 1
37 S. bracteolactus 13 1 P 246.8 SS plot 1
38 Grass E 13 1 P 81.3 GE plot 1
39 Grass E 13 1 P 82.4 GE plot 1
40 Grass E 13 1 P 78.1 GE plot 1
68 M. spinosum 26 1 P 194.4 SM plot 1
69 M. spinosum 26 1 P 166.0 SM plot 1
70 M. spinosum 26 1 P 212.0 SM plot 1
71 S. bracteolactus 26 1 P 232.2 SS plot 1
72 S. bracteolactus 26 1 P 190.1 SS plot 1
73 S. bracteolactus 26 1 P 241.6 SS plot 1
74 Grass E 26 1 P 53.9 GE plot 1
75 Grass E 26 1 P 71.9 GE plot 1
76 Grass E 26 1 P 63.1 GE plot 1
104 M. spinosum 36 1 P 178.0 SM plot 1
105 M. spinosum 36 1 P 165.9 SM plot 1
106 M. spinosum 36 1 P 144.0 SM plot 1
107 S. bracteolactus 36 1 P 144.1 SS plot 1
108 S. bracteolactus 36 1 P 142.7 SS plot 1
109 S. bracteolactus 36 1 P 170.0 SS plot 1
110 Grass E 36 1 P 45.7 GE plot 1
111 Grass E 36 1 P 37.3 GE plot 1
112 Grass E 36 1 P 33.3 GE plot 1
140 M. spinosum 57 1 P 103.5 SM plot 1
141 M. spinosum 57 1 P 115.2 SM plot 1
142 M. spinosum 57 1 P 115.2 SM plot 1
143 S. bracteolactus 57 1 P 76.7 SS plot 1
144 S. bracteolactus 57 1 P 110.5 SS plot 1
145 S. bracteolactus 57 1 P 148.5 SS plot 1
146 Grass E 57 1 P 23.0 GE plot 1
147 Grass E 57 1 P 23.3 GE plot 1
148 Grass E 57 1 P 25.6 GE plot 1
176 M. spinosum 68 1 P 65.8 SM plot 1
177 M. spinosum 68 1 P 86.9 SM plot 1
178 M. spinosum 68 1 P 79.5 SM plot 1
179 S. bracteolactus 68 1 P 70.3 SS plot 1
180 S. bracteolactus 68 1 P 67.5 SS plot 1
181 S. bracteolactus 68 1 P 87.6 SS plot 1
182 Grass E 68 1 P 14.0 GE plot 1
183 Grass E 68 1 P 24.3 GE plot 1
184 Grass E 68 1 P 17.1 GE plot 1
212 M. spinosum 80 1 P 67.0 SM plot 1
213 M. spinosum 80 1 P 70.5 SM plot 1
214 M. spinosum 80 1 P 57.9 SM plot 1
215 S. bracteolactus 80 1 P 70.6 SS plot 1
216 S. bracteolactus 80 1 P 55.4 SS plot 1
217 S. bracteolactus 80 1 P 86.0 SS plot 1
218 Grass E 80 1 P 11.4 GE plot 1
219 Grass E 80 1 P 16.2 GE plot 1
220 Grass E 80 1 P 12.8 GE plot 1
9 M. spinosum 1 2 P 251.7 SM plot 2
10 M. spinosum 1 2 P 238.7 SM plot 2
11 M. spinosum 1 2 P 252.6 SM plot 2
12 S. bracteolactus 1 2 P 257.6 SS plot 2
13 S. bracteolactus 1 2 P 265.7 SS plot 2
14 S. bracteolactus 1 2 P 273.9 SS plot 2
15 Grass E 1 2 P 60.6 GE plot 2
16 Grass E 1 2 P 53.7 GE plot 2
41 M. spinosum 13 2 P 228.0 SM plot 2
42 M. spinosum 13 2 P 231.4 SM plot 2
43 M. spinosum 13 2 P 226.6 SM plot 2
44 S. bracteolactus 13 2 P 222.7 SS plot 2
45 S. bracteolactus 13 2 P 243.9 SS plot 2
46 S. bracteolactus 13 2 P 218.6 SS plot 2
47 Grass E 13 2 P 50.5 GE plot 2
48 Grass E 13 2 P 64.3 GE plot 2
49 Grass E 13 2 P 64.3 GE plot 2
77 M. spinosum 26 2 P 185.8 SM plot 2
78 M. spinosum 26 2 P 169.7 SM plot 2
79 M. spinosum 26 2 P 154.8 SM plot 2
80 S. bracteolactus 26 2 P 207.2 SS plot 2
81 S. bracteolactus 26 2 P 200.0 SS plot 2
82 S. bracteolactus 26 2 P 204.9 SS plot 2
83 Grass E 26 2 P 50.5 GE plot 2
84 Grass E 26 2 P 43.2 GE plot 2
85 Grass E 26 2 P 52.7 GE plot 2
113 M. spinosum 36 2 P 115.9 SM plot 2
114 M. spinosum 36 2 P 124.4 SM plot 2
115 M. spinosum 36 2 P 130.7 SM plot 2
116 S. bracteolactus 36 2 P 214.9 SS plot 2
117 S. bracteolactus 36 2 P 159.2 SS plot 2
118 S. bracteolactus 36 2 P 135.9 SS plot 2
119 Grass E 36 2 P 32.1 GE plot 2
120 Grass E 36 2 P 34.0 GE plot 2
121 Grass E 36 2 P 39.2 GE plot 2
149 M. spinosum 57 2 P 96.0 SM plot 2
150 M. spinosum 57 2 P 76.7 SM plot 2
151 M. spinosum 57 2 P 101.2 SM plot 2
152 S. bracteolactus 57 2 P 80.6 SS plot 2
153 S. bracteolactus 57 2 P 83.7 SS plot 2
154 S. bracteolactus 57 2 P 108.0 SS plot 2
155 Grass E 57 2 P 20.4 GE plot 2
156 Grass E 57 2 P 23.7 GE plot 2
157 Grass E 57 2 P 21.7 GE plot 2
185 M. spinosum 68 2 P 62.2 SM plot 2
186 M. spinosum 68 2 P 70.6 SM plot 2
187 M. spinosum 68 2 P 76.6 SM plot 2
188 S. bracteolactus 68 2 P 72.8 SS plot 2
189 S. bracteolactus 68 2 P 85.2 SS plot 2
190 S. bracteolactus 68 2 P 75.4 SS plot 2
191 Grass E 68 2 P 15.4 GE plot 2
192 Grass E 68 2 P 11.8 GE plot 2
193 Grass E 68 2 P 15.3 GE plot 2
221 M. spinosum 80 2 P 38.1 SM plot 2
222 M. spinosum 80 2 P 61.2 SM plot 2
223 M. spinosum 80 2 P 53.2 SM plot 2
224 S. bracteolactus 80 2 P 72.2 SS plot 2
225 S. bracteolactus 80 2 P 64.2 SS plot 2
226 S. bracteolactus 80 2 P 61.4 SS plot 2
227 Grass E 80 2 P 13.9 GE plot 2
228 Grass E 80 2 P 11.2 GE plot 2
229 Grass E 80 2 P 13.3 GE plot 2
17 M. spinosum 1 3 P 224.3 SM plot 3
18 M. spinosum 1 3 P 207.3 SM plot 3
19 S. bracteolactus 1 3 P 240.3 SS plot 3
20 S. bracteolactus 1 3 P 215.8 SS plot 3
21 Grass E 1 3 P 75.3 GE plot 3
22 Grass E 1 3 P 65.7 GE plot 3
50 M. spinosum 13 3 P 240.8 SM plot 3
51 M. spinosum 13 3 P 178.4 SM plot 3
52 M. spinosum 13 3 P 231.2 SM plot 3
53 S. bracteolactus 13 3 P 218.6 SS plot 3
54 S. bracteolactus 13 3 P 244.1 SS plot 3
55 S. bracteolactus 13 3 P 204.4 SS plot 3
56 Grass E 13 3 P 69.3 GE plot 3
57 Grass E 13 3 P 77.2 GE plot 3
58 Grass E 13 3 P 62.5 GE plot 3
86 M. spinosum 26 3 P 154.7 SM plot 3
87 M. spinosum 26 3 P 221.0 SM plot 3
88 M. spinosum 26 3 P 231.7 SM plot 3
89 S. bracteolactus 26 3 P 171.4 SS plot 3
90 S. bracteolactus 26 3 P 209.4 SS plot 3
91 S. bracteolactus 26 3 P 180.5 SS plot 3
92 Grass E 26 3 P 54.7 GE plot 3
93 Grass E 26 3 P 45.9 GE plot 3
94 Grass E 26 3 P 41.1 GE plot 3
122 M. spinosum 36 3 P 164.7 SM plot 3
123 M. spinosum 36 3 P 142.9 SM plot 3
124 M. spinosum 36 3 P 161.4 SM plot 3
125 S. bracteolactus 36 3 P 152.0 SS plot 3
126 S. bracteolactus 36 3 P 114.6 SS plot 3
127 S. bracteolactus 36 3 P 183.8 SS plot 3
128 Grass E 36 3 P 24.4 GE plot 3
129 Grass E 36 3 P 28.0 GE plot 3
130 Grass E 36 3 P 26.4 GE plot 3
158 M. spinosum 57 3 P 129.1 SM plot 3
159 M. spinosum 57 3 P 119.9 SM plot 3
160 M. spinosum 57 3 P 148.1 SM plot 3
161 S. bracteolactus 57 3 P 110.6 SS plot 3
162 S. bracteolactus 57 3 P 54.1 SS plot 3
163 S. bracteolactus 57 3 P 95.1 SS plot 3
164 Grass E 57 3 P 31.3 GE plot 3
165 Grass E 57 3 P 22.5 GE plot 3
166 Grass E 57 3 P 27.4 GE plot 3
[ reached 'max' / getOption("max.print") -- omitted 81 rows ]
lfmc.gd$plot <- with(lfmc.gd, factor(plot, levels = c("4", "5", "6", "1", "2", "3")))
lfmc.gd
Grouped Data: lfmc ~ time | group
leaf.type time plot site lfmc group
1 M. spinosum 1 1 P 238.4 SM plot 1
2 M. spinosum 1 1 P 269.8 SM plot 1
3 M. spinosum 1 1 P 325.2 SM plot 1
4 S. bracteolactus 1 1 P 326.7 SS plot 1
5 S. bracteolactus 1 1 P 257.3 SS plot 1
6 S. bracteolactus 1 1 P 279.3 SS plot 1
7 Grass E 1 1 P 91.9 GE plot 1
8 Grass E 1 1 P 95.0 GE plot 1
32 M. spinosum 13 1 P 215.8 SM plot 1
33 M. spinosum 13 1 P 216.7 SM plot 1
34 M. spinosum 13 1 P 205.3 SM plot 1
35 S. bracteolactus 13 1 P 208.5 SS plot 1
36 S. bracteolactus 13 1 P 225.7 SS plot 1
37 S. bracteolactus 13 1 P 246.8 SS plot 1
38 Grass E 13 1 P 81.3 GE plot 1
39 Grass E 13 1 P 82.4 GE plot 1
40 Grass E 13 1 P 78.1 GE plot 1
68 M. spinosum 26 1 P 194.4 SM plot 1
69 M. spinosum 26 1 P 166.0 SM plot 1
70 M. spinosum 26 1 P 212.0 SM plot 1
71 S. bracteolactus 26 1 P 232.2 SS plot 1
72 S. bracteolactus 26 1 P 190.1 SS plot 1
73 S. bracteolactus 26 1 P 241.6 SS plot 1
74 Grass E 26 1 P 53.9 GE plot 1
75 Grass E 26 1 P 71.9 GE plot 1
76 Grass E 26 1 P 63.1 GE plot 1
104 M. spinosum 36 1 P 178.0 SM plot 1
105 M. spinosum 36 1 P 165.9 SM plot 1
106 M. spinosum 36 1 P 144.0 SM plot 1
107 S. bracteolactus 36 1 P 144.1 SS plot 1
108 S. bracteolactus 36 1 P 142.7 SS plot 1
109 S. bracteolactus 36 1 P 170.0 SS plot 1
110 Grass E 36 1 P 45.7 GE plot 1
111 Grass E 36 1 P 37.3 GE plot 1
112 Grass E 36 1 P 33.3 GE plot 1
140 M. spinosum 57 1 P 103.5 SM plot 1
141 M. spinosum 57 1 P 115.2 SM plot 1
142 M. spinosum 57 1 P 115.2 SM plot 1
143 S. bracteolactus 57 1 P 76.7 SS plot 1
144 S. bracteolactus 57 1 P 110.5 SS plot 1
145 S. bracteolactus 57 1 P 148.5 SS plot 1
146 Grass E 57 1 P 23.0 GE plot 1
147 Grass E 57 1 P 23.3 GE plot 1
148 Grass E 57 1 P 25.6 GE plot 1
176 M. spinosum 68 1 P 65.8 SM plot 1
177 M. spinosum 68 1 P 86.9 SM plot 1
178 M. spinosum 68 1 P 79.5 SM plot 1
179 S. bracteolactus 68 1 P 70.3 SS plot 1
180 S. bracteolactus 68 1 P 67.5 SS plot 1
181 S. bracteolactus 68 1 P 87.6 SS plot 1
182 Grass E 68 1 P 14.0 GE plot 1
183 Grass E 68 1 P 24.3 GE plot 1
184 Grass E 68 1 P 17.1 GE plot 1
212 M. spinosum 80 1 P 67.0 SM plot 1
213 M. spinosum 80 1 P 70.5 SM plot 1
214 M. spinosum 80 1 P 57.9 SM plot 1
215 S. bracteolactus 80 1 P 70.6 SS plot 1
216 S. bracteolactus 80 1 P 55.4 SS plot 1
217 S. bracteolactus 80 1 P 86.0 SS plot 1
218 Grass E 80 1 P 11.4 GE plot 1
219 Grass E 80 1 P 16.2 GE plot 1
220 Grass E 80 1 P 12.8 GE plot 1
9 M. spinosum 1 2 P 251.7 SM plot 2
10 M. spinosum 1 2 P 238.7 SM plot 2
11 M. spinosum 1 2 P 252.6 SM plot 2
12 S. bracteolactus 1 2 P 257.6 SS plot 2
13 S. bracteolactus 1 2 P 265.7 SS plot 2
14 S. bracteolactus 1 2 P 273.9 SS plot 2
15 Grass E 1 2 P 60.6 GE plot 2
16 Grass E 1 2 P 53.7 GE plot 2
41 M. spinosum 13 2 P 228.0 SM plot 2
42 M. spinosum 13 2 P 231.4 SM plot 2
43 M. spinosum 13 2 P 226.6 SM plot 2
44 S. bracteolactus 13 2 P 222.7 SS plot 2
45 S. bracteolactus 13 2 P 243.9 SS plot 2
46 S. bracteolactus 13 2 P 218.6 SS plot 2
47 Grass E 13 2 P 50.5 GE plot 2
48 Grass E 13 2 P 64.3 GE plot 2
49 Grass E 13 2 P 64.3 GE plot 2
77 M. spinosum 26 2 P 185.8 SM plot 2
78 M. spinosum 26 2 P 169.7 SM plot 2
79 M. spinosum 26 2 P 154.8 SM plot 2
80 S. bracteolactus 26 2 P 207.2 SS plot 2
81 S. bracteolactus 26 2 P 200.0 SS plot 2
82 S. bracteolactus 26 2 P 204.9 SS plot 2
83 Grass E 26 2 P 50.5 GE plot 2
84 Grass E 26 2 P 43.2 GE plot 2
85 Grass E 26 2 P 52.7 GE plot 2
113 M. spinosum 36 2 P 115.9 SM plot 2
114 M. spinosum 36 2 P 124.4 SM plot 2
115 M. spinosum 36 2 P 130.7 SM plot 2
116 S. bracteolactus 36 2 P 214.9 SS plot 2
117 S. bracteolactus 36 2 P 159.2 SS plot 2
118 S. bracteolactus 36 2 P 135.9 SS plot 2
119 Grass E 36 2 P 32.1 GE plot 2
120 Grass E 36 2 P 34.0 GE plot 2
121 Grass E 36 2 P 39.2 GE plot 2
149 M. spinosum 57 2 P 96.0 SM plot 2
150 M. spinosum 57 2 P 76.7 SM plot 2
151 M. spinosum 57 2 P 101.2 SM plot 2
152 S. bracteolactus 57 2 P 80.6 SS plot 2
153 S. bracteolactus 57 2 P 83.7 SS plot 2
154 S. bracteolactus 57 2 P 108.0 SS plot 2
155 Grass E 57 2 P 20.4 GE plot 2
156 Grass E 57 2 P 23.7 GE plot 2
157 Grass E 57 2 P 21.7 GE plot 2
185 M. spinosum 68 2 P 62.2 SM plot 2
186 M. spinosum 68 2 P 70.6 SM plot 2
187 M. spinosum 68 2 P 76.6 SM plot 2
188 S. bracteolactus 68 2 P 72.8 SS plot 2
189 S. bracteolactus 68 2 P 85.2 SS plot 2
190 S. bracteolactus 68 2 P 75.4 SS plot 2
191 Grass E 68 2 P 15.4 GE plot 2
192 Grass E 68 2 P 11.8 GE plot 2
193 Grass E 68 2 P 15.3 GE plot 2
221 M. spinosum 80 2 P 38.1 SM plot 2
222 M. spinosum 80 2 P 61.2 SM plot 2
223 M. spinosum 80 2 P 53.2 SM plot 2
224 S. bracteolactus 80 2 P 72.2 SS plot 2
225 S. bracteolactus 80 2 P 64.2 SS plot 2
226 S. bracteolactus 80 2 P 61.4 SS plot 2
227 Grass E 80 2 P 13.9 GE plot 2
228 Grass E 80 2 P 11.2 GE plot 2
229 Grass E 80 2 P 13.3 GE plot 2
17 M. spinosum 1 3 P 224.3 SM plot 3
18 M. spinosum 1 3 P 207.3 SM plot 3
19 S. bracteolactus 1 3 P 240.3 SS plot 3
20 S. bracteolactus 1 3 P 215.8 SS plot 3
21 Grass E 1 3 P 75.3 GE plot 3
22 Grass E 1 3 P 65.7 GE plot 3
50 M. spinosum 13 3 P 240.8 SM plot 3
51 M. spinosum 13 3 P 178.4 SM plot 3
52 M. spinosum 13 3 P 231.2 SM plot 3
53 S. bracteolactus 13 3 P 218.6 SS plot 3
54 S. bracteolactus 13 3 P 244.1 SS plot 3
55 S. bracteolactus 13 3 P 204.4 SS plot 3
56 Grass E 13 3 P 69.3 GE plot 3
57 Grass E 13 3 P 77.2 GE plot 3
58 Grass E 13 3 P 62.5 GE plot 3
86 M. spinosum 26 3 P 154.7 SM plot 3
87 M. spinosum 26 3 P 221.0 SM plot 3
88 M. spinosum 26 3 P 231.7 SM plot 3
89 S. bracteolactus 26 3 P 171.4 SS plot 3
90 S. bracteolactus 26 3 P 209.4 SS plot 3
91 S. bracteolactus 26 3 P 180.5 SS plot 3
92 Grass E 26 3 P 54.7 GE plot 3
93 Grass E 26 3 P 45.9 GE plot 3
94 Grass E 26 3 P 41.1 GE plot 3
122 M. spinosum 36 3 P 164.7 SM plot 3
123 M. spinosum 36 3 P 142.9 SM plot 3
124 M. spinosum 36 3 P 161.4 SM plot 3
125 S. bracteolactus 36 3 P 152.0 SS plot 3
126 S. bracteolactus 36 3 P 114.6 SS plot 3
127 S. bracteolactus 36 3 P 183.8 SS plot 3
128 Grass E 36 3 P 24.4 GE plot 3
129 Grass E 36 3 P 28.0 GE plot 3
130 Grass E 36 3 P 26.4 GE plot 3
158 M. spinosum 57 3 P 129.1 SM plot 3
159 M. spinosum 57 3 P 119.9 SM plot 3
160 M. spinosum 57 3 P 148.1 SM plot 3
161 S. bracteolactus 57 3 P 110.6 SS plot 3
162 S. bracteolactus 57 3 P 54.1 SS plot 3
163 S. bracteolactus 57 3 P 95.1 SS plot 3
164 Grass E 57 3 P 31.3 GE plot 3
165 Grass E 57 3 P 22.5 GE plot 3
166 Grass E 57 3 P 27.4 GE plot 3
[ reached 'max' / getOption("max.print") -- omitted 81 rows ]
ggplot(lfmc, aes(time, lfmc)) +
geom_point() +
geom_smooth(method = "loess") +
facet_wrap(~leaf.type) +
xlab("Time (days)") +
ylab("LMFC (%)")
`geom_smooth()` using formula 'y ~ x'
nlsL <- nlsList(lfmc ~ SSdlf(time, A, w, m, s), data = lfmc.gd)
Warning: 3 errors caught in nls(model, data = data, control = controlvals). The error messages and their frequencies are
step factor 0.000488281 reduced below 'minFactor' of 0.000976562
1
number of iterations exceeded maximum of 50
2
coef(nlsL)
plot(intervals(nlsL))
nl.re <- nlme(nlsL, control = nlmeControl(maxIter = 5000, msMaxIter = 1500))
Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
Iteration 1, LME step: nlminb() did not converge (code = 1). PORT message: function evaluation limit reached without convergence (9)
nl.re.1 <- update(nl.re, random = pdDiag(A + w + m + s ~ 1))
fxf <- fixef(nl.re.1)
fxf
A w m s
193.45267 21.45740 31.48270 -22.59711
# nl.me <- update(nl.re.1, fixed = list(A + w + m + s ~ leaf.type),
# start = c(A = fxf[1],0,0,0,
# w = fxf[2],0,0,0,
# m = fxf[3],0,0,0,
# s = fxf[4],0,0,0))
nl.re.2 <- update(nl.re.1, random = pdDiag(A + w + m ~ 1))
nl.re.2
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -1070.648
Fixed: list(A ~ 1, w ~ 1, m ~ 1, s ~ 1)
A w m s
193.46656 21.43199 31.48851 -22.60602
Random effects:
Formula: list(A ~ 1, w ~ 1, m ~ 1)
Level: group
Structure: Diagonal
A w m Residual
StdDev: 119.3792 10.94166 8.927276 15.26574
Number of Observations: 247
Number of Groups: 12
fxf <- fixef(nl.re.2)
fxf
A w m s
193.46656 21.43199 31.48851 -22.60602
nl.me.1 <- update(
nl.re.2,
fixed = list(A + w + m + s ~ leaf.type),
start = c(A = fxf[1], 0, 0, 0, w = fxf[2], 0, 0, 0, m = fxf[3], 0, 0, 0, s = fxf[4], 0, 0, 0)
)
nl.me.1
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -1034.375
Fixed: list(A + w + m + s ~ leaf.type)
A.(Intercept) A.leaf.typeGrass E A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
23.135414 53.294908 296.620596 263.313118
w.(Intercept) w.leaf.typeGrass E w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
50.291577 -34.313924 -26.732287 2.135665
m.(Intercept) m.leaf.typeGrass E m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus
51.533983 -21.920404 -20.647403 -18.170359
s.(Intercept) s.leaf.typeGrass E s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus
17.605794 -25.926031 -43.529450 -33.666226
Random effects:
Formula: list(A ~ 1, w ~ 1, m ~ 1)
Level: group
Structure: Diagonal
A.(Intercept) w.(Intercept) m.(Intercept) Residual
StdDev: 8.112813 0.001236897 2.468334 15.38607
Number of Observations: 247
Number of Groups: 12
AIC(nl.re.2, nl.me.1)
IC_tab(nl.re.2, nl.me.1)
The small StdDev of “w” suggests to remove the random-effect of plot on this parameter.
nl.me.2 <- update(nl.me.1, random = pdDiag(A + m ~ 1))
IC_tab(nl.me.1, nl.me.2)
in addition, if the random-effects on “m” is also removed.
nl.me.3 <- update(nl.me.2, random = pdDiag(A ~ 1))
nl.me.3
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -1035.19
Fixed: list(A + w + m + s ~ leaf.type)
A.(Intercept) A.leaf.typeGrass E A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
23.071500 53.950821 408.256907 258.698322
w.(Intercept) w.leaf.typeGrass E w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
50.311106 -34.625804 -79.299174 3.717378
m.(Intercept) m.leaf.typeGrass E m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus
51.608221 -21.975396 -32.109785 -17.665025
s.(Intercept) s.leaf.typeGrass E s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus
17.683236 -26.518320 -60.147105 -33.113726
Random effects:
Formula: A ~ 1 | group
A.(Intercept) Residual
StdDev: 10.6035 15.53457
Number of Observations: 247
Number of Groups: 12
the model fit is similar, supporting the simpler model.
IC_tab(nl.me.2, nl.me.3)
plot(nl.me.3)
hist(residuals(nl.me.3, type = "normalized"), breaks = 20)
rse <- residuals(nl.me.3) / sigma(nl.me.3)
par(pty = "s")
qqnorm(rse)
qqline(rse)
Heterocedasticity: variance modeling —-
nl.me.3.vm <- update(nl.me.3, weights = varIdent(form = ~ 1 | leaf.type))
nl.me.3.vm
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -965.5628
Fixed: list(A + w + m + s ~ leaf.type)
A.(Intercept) A.leaf.typeGrass E A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
50.1522851 28.9553861 346.5816334 231.1569667
w.(Intercept) w.leaf.typeGrass E w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
24.6693607 -9.4837876 -35.8327942 29.7533014
m.(Intercept) m.leaf.typeGrass E m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus
48.9613646 -19.9780200 -26.2920794 -14.9909275
s.(Intercept) s.leaf.typeGrass E s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus
-16.2294127 6.3520529 -21.0664455 0.9389323
Random effects:
Formula: A ~ 1 | group
A.(Intercept) Residual
StdDev: 9.412263 6.824811
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.0000000 0.8547835 3.1456077 3.0262116
Number of Observations: 247
Number of Groups: 12
IC_tab(nl.me.3, nl.me.3.vm)
여기서 가중잔차는 weight를 이용해서 계산해야 제대로 plot이 됨
res <-
tibble(
fitted = fitted(nl.me.3.vm),
resid = resid(nl.me.3.vm), # 일반 잔차
nresid = resid(nl.me.3.vm, type = "normalized"), # 표준 잔차
weight = diag(var_cov(nl.me.3.vm)), # 가중치
wresid = sqrt(1 / weight) * resid # 가중잔차 수동계산
) %>%
bind_cols(lfmc.gd)
res
res %>%
ggplot(aes(fitted, nresid)) +
geom_point(aes(color = leaf.type))
res %>%
ggplot(aes(fitted, wresid)) +
geom_point(aes(color = leaf.type))
plot(nl.me.3.vm)
hist(residuals(nl.me.3.vm, type = "normalized"), main = "", xlab = "Residuals")
the residuals look fine now.
qqnorm(resid(nl.me.3.vm, type = "normalized"))
qqline(resid(nl.me.3.vm, type = "normalized"))
Temporal correlation check. This plot suggests there is no need for modeling the correlation of the residuals.
plot(ACF(nl.me.3.vm), alpha = 0.01)
nl.me.ml <- update(nl.me.3.vm, method = "ML")
nl.me.ml
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -965.5628
Fixed: list(A + w + m + s ~ leaf.type)
A.(Intercept) A.leaf.typeGrass E A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
50.1522851 28.9553861 346.5816334 231.1569667
w.(Intercept) w.leaf.typeGrass E w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
24.6693607 -9.4837876 -35.8327942 29.7533014
m.(Intercept) m.leaf.typeGrass E m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus
48.9613646 -19.9780200 -26.2920794 -14.9909275
s.(Intercept) s.leaf.typeGrass E s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus
-16.2294127 6.3520529 -21.0664455 0.9389323
Random effects:
Formula: A ~ 1 | group
A.(Intercept) Residual
StdDev: 9.412263 6.824811
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.0000000 0.8547835 3.1456077 3.0262116
Number of Observations: 247
Number of Groups: 12
the model is reduced removing the fixed-effect of “leaf type” on “s”
nl.me.ml.1 <- update(nl.me.ml,
fixed = list(A + w + m ~ leaf.type, s ~ 1),
start = c(A = fxf[1], 0, 0, 0, w = fxf[2], 0, 0, 0, m = fxf[3], 0, 0, 0, s = fxf[4])
)
nl.me.ml.1
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -969.2617
Fixed: list(A + w + m ~ leaf.type, s ~ 1)
A.(Intercept) A.leaf.typeGrass E A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
49.94562 39.60979 219.45652 232.11909
w.(Intercept) w.leaf.typeGrass E w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
25.24501 -14.34293 35.36122 28.70389
m.(Intercept) m.leaf.typeGrass E m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus
48.30390 -21.59477 -15.57765 -14.37454
s
-15.44291
Random effects:
Formula: A ~ 1 | group
A.(Intercept) Residual
StdDev: 9.612472 6.825604
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.0000000 0.8714808 3.2750696 3.0238816
Number of Observations: 247
Number of Groups: 12
the AIC comparison suggests that the more complex model is not justified (fits are similar)
IC_tab(nl.me.ml, nl.me.ml.1)
the model is reduced removing the fixed-effect of “leaf type” on “m”
nl.me.ml.2 <- update(nl.me.ml,
fixed = list(A + w ~ leaf.type, m + s ~ 1),
start = c(A = fxf[1], 0, 0, 0, w = fxf[2], 0, 0, 0, m = fxf[3], s = fxf[4])
)
nl.me.ml.2
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -972.7465
Fixed: list(A + w ~ leaf.type, m + s ~ 1)
A.(Intercept) A.leaf.typeGrass E A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
54.25730 30.93396 223.27914 240.31542
w.(Intercept) w.leaf.typeGrass E w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
29.05803 -20.68663 31.68646 26.99003
m s
30.91655 -16.15435
Random effects:
Formula: A ~ 1 | group
A.(Intercept) Residual
StdDev: 9.216113 7.028345
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.0000000 0.8820414 3.1693973 2.9643303
Number of Observations: 247
Number of Groups: 12
the model is reduced removing the fixed-effect of “leaf type” on “W”.
nl.me.ml.3 <- update(nl.me.ml,
fixed = list(A ~ leaf.type, w + m + s ~ 1),
start = c(A = fxf[1], 0, 0, 0, w = fxf[2], m = fxf[3], s = fxf[4])
)
nl.me.ml.3
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -1017.91
Fixed: list(A ~ leaf.type, w + m + s ~ 1)
A.(Intercept) A.leaf.typeGrass E A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
63.92200 14.27252 251.99435 265.42764
w m s
17.33019 32.65673 -25.16298
Random effects:
Formula: A ~ 1 | group
A.(Intercept) Residual
StdDev: 8.138947 7.795662
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.000000 1.584887 2.813985 2.731415
Number of Observations: 247
Number of Groups: 12
nl.me.ml.4 <- update(nl.me.ml,
fixed = list(w ~ leaf.type, A + m + s ~ 1),
start = c(A = fxf[1], w = fxf[2], 0, 0, 0, m = fxf[3], s = fxf[4])
)
nl.me.ml.4
Nonlinear mixed-effects model fit by maximum likelihood
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
Log-likelihood: -999.8224
Fixed: list(w ~ leaf.type, A + m + s ~ 1)
w.(Intercept) w.leaf.typeGrass E w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
29.16029 -20.38028 33.19430 28.72481
A m s
175.78131 31.08598 -15.70297
Random effects:
Formula: A ~ 1 | group
A Residual
StdDev: 108.0655 7.042308
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.0000000 0.8770721 3.1889180 2.9304615
Number of Observations: 247
Number of Groups: 12
IC_tab(nl.me.ml.1, nl.me.ml.2, nl.me.ml.3, nl.me.ml.4)
This model assumes:
M1 <- update(nl.me.ml.2, method = "REML")
summary(M1)
Nonlinear mixed-effects model fit by REML
Model: lfmc ~ SSdlf(time, A, w, m, s)
Data: lfmc.gd
AIC BIC logLik
1931.668 1983.689 -950.834
Random effects:
Formula: A ~ 1 | group
A.(Intercept) Residual
StdDev: 9.408521 7.162899
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W Grass E M. spinosum S. bracteolactus
1.0000000 0.8851764 3.1796102 2.9647393
Fixed effects: list(A + w ~ leaf.type, m + s ~ 1)
Value Std.Error DF t-value p-value
A.(Intercept) 54.25350 5.988181 226 9.060097 0e+00
A.leaf.typeGrass E 30.92293 8.835210 226 3.499966 6e-04
A.leaf.typeM. spinosum 223.24504 15.915557 226 14.026844 0e+00
A.leaf.typeS. bracteolactus 240.27654 16.857356 226 14.253513 0e+00
w.(Intercept) 29.05581 1.712960 226 16.962334 0e+00
w.leaf.typeGrass E -20.68938 2.592469 226 -7.980569 0e+00
w.leaf.typeM. spinosum 31.67297 7.753176 226 4.085161 1e-04
w.leaf.typeS. bracteolactus 26.97508 8.071111 226 3.342177 1e-03
m 30.92881 2.065196 226 14.976213 0e+00
s -16.15406 2.365392 226 -6.829339 0e+00
Correlation:
A.(In) A.l.GE A..M.s A..S.b w.(In) w.l.GE w..M.s w..S.b m
A.leaf.typeGrass E -0.527
A.leaf.typeM. spinosum -0.146 0.535
A.leaf.typeS. bracteolactus -0.115 0.537 0.746
w.(Intercept) -0.217 -0.034 -0.200 -0.216
w.leaf.typeGrass E -0.035 -0.283 -0.382 -0.399 -0.258
w.leaf.typeM. spinosum -0.118 -0.227 -0.547 -0.454 0.157 0.573
w.leaf.typeS. bracteolactus -0.129 -0.241 -0.462 -0.575 0.188 0.601 0.640
m -0.217 -0.324 -0.633 -0.666 0.092 0.145 0.161 0.172
s -0.246 -0.354 -0.710 -0.748 0.416 0.575 0.702 0.752 0.544
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.27999661 -0.59657156 -0.02042043 0.57273326 3.12018001
Number of Observations: 247
Number of Groups: 12
fixef(M1)
A.(Intercept) A.leaf.typeGrass E A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus
54.25350 30.92293 223.24504 240.27654
w.(Intercept) w.leaf.typeGrass E w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus
29.05581 -20.68938 31.67297 26.97508
m s
30.92881 -16.15406
ranef(M1)
A.(Intercept)
GW plot 4 -1.915443
GW plot 5 -1.923198
GW plot 6 3.838641
GE plot 1 15.652459
SM plot 1 6.967303
SS plot 1 9.166925
GE plot 2 -11.066197
SM plot 2 -5.358910
SS plot 2 2.442781
GE plot 3 -4.586263
SM plot 3 -1.608393
SS plot 3 -11.609707
intervals(M1)
Approximate 95% confidence intervals
Fixed effects:
lower est. upper
A.(Intercept) 42.45370 54.25350 66.05331
A.leaf.typeGrass E 13.51301 30.92293 48.33286
A.leaf.typeM. spinosum 191.88318 223.24504 254.60691
A.leaf.typeS. bracteolactus 207.05884 240.27654 273.49423
w.(Intercept) 25.68039 29.05581 32.43122
w.leaf.typeGrass E -25.79788 -20.68938 -15.58088
w.leaf.typeM. spinosum 16.39521 31.67297 46.95073
w.leaf.typeS. bracteolactus 11.07083 26.97508 42.87934
m 26.85931 30.92881 34.99831
s -20.81510 -16.15406 -11.49302
Random Effects:
Level: group
lower est. upper
sd(A.(Intercept)) 5.0352 9.408521 17.58029
Variance function:
lower est. upper
Grass E 0.675418 0.8851764 1.160077
M. spinosum 2.451780 3.1796102 4.123502
S. bracteolactus 2.285695 2.9647393 3.845517
Within-group standard error:
lower est. upper
5.966409 7.162899 8.599330
# 2.2.1 - Overall predictions (Table 3) ----
Agw <- fixef(M1)[1] # "A" for grasses in the W site (GW)
Age <- fixef(M1)[1] + fixef(M1)[2] # "A" for grasses in the E site (GE)
Asm <- fixef(M1)[1] + fixef(M1)[3] # "A" for Mullinum spinosum (SM)
Ass <- fixef(M1)[1] + fixef(M1)[4] # "A" for Senecio filaginoides (SS)
Wgw <- fixef(M1)[5] # "w" for grasses in the W site (GW)
Wge <- fixef(M1)[5] + fixef(M1)[6] # "w" for grasses in the E site (GE)
Wsm <- fixef(M1)[5] + fixef(M1)[7] # "w" for Mullinum spinosum (SM)
Wss <- fixef(M1)[5] + fixef(M1)[8] # "w" for Senecio filaginoides (SS)
M <- fixef(M1)[9] # "m" for all the leaf types
S <- fixef(M1)[10] # "s" for all the leaf types
GW <- subset(lfmc, leaf.type == "Grass W")
GE <- subset(lfmc, leaf.type == "Grass E")
SM <- subset(lfmc, leaf.type == "M. spinosum")
SS <- subset(lfmc, leaf.type == "S. bracteolactus")
par(mfrow = c(2, 2), cex = 0.75)
# Grasses W
plot(GW$time, GW$lfmc, xlab = "", ylab = "LFMC (%)", ylim = c(0, 100), main = "Grasses W site", col = "darkgreen")
curve(((Agw - Wgw) / (1 + exp((M - x) / S)) + Wgw), lwd = 2, lty = 1, add = T, col = "darkgreen")
# Grasses E
plot(GE$time, GE$lfmc, xlab = "", ylab = "LFMC (%)", ylim = c(0, 100), main = "Grasses E site", col = "green")
curve(((Age - Wge) / (1 + exp((M - x) / S)) + Wge), lwd = 2, lty = 1, add = T, col = "green")
# M. Spinosum
plot(SM$time, SM$lfmc, xlab = "Time (days)", ylab = "LFMC (%)", ylim = c(0, 350), main = "M. spinosum", font.main = 4, col = "darkorange")
curve(((Asm - Wsm) / (1 + exp((M - x) / S)) + Wsm), lwd = 2, lty = 1, add = T, col = "darkorange")
# S. filaginoides
plot(SS$time, SS$lfmc, xlab = "Time (days)", ylab = "LFMC (%)", ylim = c(0, 350), main = "S. bracteolactus", font.main = 4, col = "darkred")
curve(((Ass - Wss) / (1 + exp((M - x) / S)) + Wss), lwd = 2, lty = 1, add = T, col = "darkred")
ranef(M1)
A.(Intercept)
GW plot 4 -1.915443
GW plot 5 -1.923198
GW plot 6 3.838641
GE plot 1 15.652459
SM plot 1 6.967303
SS plot 1 9.166925
GE plot 2 -11.066197
SM plot 2 -5.358910
SS plot 2 2.442781
GE plot 3 -4.586263
SM plot 3 -1.608393
SS plot 3 -11.609707
# 2.2.2 - Predictions at plot level (Table 3) ----
Agw.p4 <- fixef(M1)[1] + ranef(M1)[1, 1] # "A" for grasses in the plot 4 (GW)
Agw.p5 <- fixef(M1)[1] + ranef(M1)[2, 1] # "A" for grasses in the plot 5 (GW)
Agw.p6 <- fixef(M1)[1] + ranef(M1)[3, 1] # "A" for grasses in the plot 6 (GW)
Age.p1 <- fixef(M1)[1] + fixef(M1)[2] + ranef(M1)[4, 1] # "A" for grasses in the plot 1 (GW)
Age.p2 <- fixef(M1)[1] + fixef(M1)[2] + ranef(M1)[7, 1] # "A" for grasses in the plot 2 (GW)
Age.p3 <- fixef(M1)[1] + fixef(M1)[2] + ranef(M1)[10, 1] # "A" for grasses in the plot 3 (GW)
Asm.p1 <- fixef(M1)[1] + fixef(M1)[3] + ranef(M1)[5, 1] # "A" for M. spinosum in the plot 1 (SM)
Asm.p2 <- fixef(M1)[1] + fixef(M1)[3] + ranef(M1)[8, 1] # "A" for M. spinosum in the plot 2 (SM)
Asm.p3 <- fixef(M1)[1] + fixef(M1)[3] + ranef(M1)[11, 1] # "A" for M. spinosum in the plot 3 (SM)
Ass.p1 <- fixef(M1)[1] + fixef(M1)[4] + ranef(M1)[6, 1] # "A" for S. filaginoides in the plot 1 (SM)
Ass.p2 <- fixef(M1)[1] + fixef(M1)[4] + ranef(M1)[9, 1] # "A" for S. filaginoides in the plot 2 (SM)
Ass.p3 <- fixef(M1)[1] + fixef(M1)[4] + ranef(M1)[12, 1] # "A" for S. filaginoides in the plot 3 (SM)
par(mfrow = c(2, 2), cex = 0.75) # (Figure 4)
# Grasses W site
plot(GW$time, GW$lfmc, xlab = "", ylab = "LFMC (%)", ylim = c(0, 100), main = "Grasses W site", col = "darkgreen")
curve(((Agw.p4 - Wgw) / (1 + exp((M - x) / S)) + Wgw), lwd = 1, lty = 2, add = T, col = "darkgreen") # plot 4
curve(((Agw.p5 - Wgw) / (1 + exp((M - x) / S)) + Wgw), lwd = 1, lty = 2, add = T, col = "darkgreen") # plot 5
curve(((Agw.p6 - Wgw) / (1 + exp((M - x) / S)) + Wgw), lwd = 1, lty = 2, add = T, col = "darkgreen") # plot 6
# Grasses E site
plot(GE$time, GE$lfmc, xlab = "", ylab = "LFMC (%)", ylim = c(0, 100), main = "Grasses E site", col = "green")
curve(((Age.p1 - Wge) / (1 + exp((M - x) / S)) + Wge), lwd = 1, lty = 2, add = T, col = "green") # plot 1
curve(((Age.p2 - Wge) / (1 + exp((M - x) / S)) + Wge), lwd = 1, lty = 2, add = T, col = "green") # plot 2
curve(((Age.p3 - Wge) / (1 + exp((M - x) / S)) + Wge), lwd = 1, lty = 2, add = T, col = "green") # plot 3
# M. Spinosum (E site)
plot(SM$time, SM$lfmc, xlab = "Time (days)", ylab = "LFMC (%)", ylim = c(0, 350), main = "M. spinosum", font.main = 4, col = "darkorange")
curve(((Asm.p1 - Wsm) / (1 + exp((M - x) / S)) + Wsm), lwd = 1, lty = 2, add = T, col = "darkorange") # plot 1
curve(((Asm.p2 - Wsm) / (1 + exp((M - x) / S)) + Wsm), lwd = 1, lty = 2, add = T, col = "darkorange") # plot 2
curve(((Asm.p3 - Wsm) / (1 + exp((M - x) / S)) + Wsm), lwd = 1, lty = 2, add = T, col = "darkorange") # plot 3
# S. filaginoides (E site)
plot(SS$time, SS$lfmc, xlab = "Time (days)", ylab = "LFMC (%)", ylim = c(0, 350), main = "S. bracteolactus", font.main = 4, col = "darkred")
curve(((Ass.p1 - Wss) / (1 + exp((M - x) / S)) + Wss), lwd = 1, lty = 2, add = T, col = "darkred") # plot 1
curve(((Ass.p2 - Wss) / (1 + exp((M - x) / S)) + Wss), lwd = 1, lty = 2, add = T, col = "darkred") # plot 2
curve(((Ass.p3 - Wss) / (1 + exp((M - x) / S)) + Wss), lwd = 1, lty = 2, add = T, col = "darkred") # plot 3
# 2.2.3 - Derivatives (drying speed) ----
lfmc.GW <- expression((Agw - Wgw) / (1 + exp((M - x) / S)) + Wgw) # LFMC dynamics for grasses in the W site
ds.GW <- deriv(lfmc.GW, "x", function.arg = TRUE) # Drying speed for grasses in the W site
lfmc.GE <- expression((Age - Wge) / (1 + exp((M - x) / S)) + Wge) # LFMC dynamics for grasses in the E site
ds.GE <- deriv(lfmc.GE, "x", function.arg = TRUE) # Drying speed for grasses in the E site
lfmc.SM <- expression((Asm - Wsm) / (1 + exp((M - x) / S)) + Wsm) # LFMC dynamics for M. spinosum
ds.SM <- deriv(lfmc.SM, "x", function.arg = TRUE) # Drying speed for M. spinosum
lfmc.SS <- expression((Ass - Wss) / (1 + exp((M - x) / S)) + Wss) # LFMC dynamics for S. filaginoides
ds.SS <- deriv(lfmc.SS, "x", function.arg = TRUE) # # Drying speed for S. filaginoides
xvec <- seq(0, 100, length = 1000)
y.ds.GW <- ds.GW(xvec)
y.ds.GE <- ds.GE(xvec)
y.ds.SM <- ds.SM(xvec)
y.ds.SS <- ds.SS(xvec)
par(mfrow = c(2, 2), cex = 0.75) # (Figure 4)
plot(xvec, y.ds.GW,
xlim = c(0, 80), ylim = c(0, 4),
ylab = "Drying speed", xlab = "Time (days)", type = "n", main = "Grasses in the W site"
)
lines(xvec, -1 * (attr(y.ds.GW, "grad")), lwd = 2, lty = 1, col = "darkgreen")
plot(xvec, y.ds.GE,
xlim = c(0, 80), ylim = c(0, 4),
ylab = "Drying speed", xlab = "Time (days)", type = "n", main = "Grasses in the E site"
)
lines(xvec, -1 * (attr(y.ds.GE, "grad")), lwd = 2, lty = 1, col = "green")
plot(xvec, y.ds.SM,
xlim = c(0, 80), ylim = c(0, 4),
ylab = "Drying speed", xlab = "Time (days)", type = "n", main = "M. spinosum", font.main = 4
)
lines(xvec, -1 * (attr(y.ds.SM, "grad")), lwd = 2, lty = 1, col = "darkorange")
plot(xvec, y.ds.SS,
xlim = c(0, 80), ylim = c(0, 4),
ylab = "Drying speed", xlab = "Time (days)", type = "n", main = "S. filaginoides", font.main = 4
)
lines(xvec, -1 * (attr(y.ds.SS, "grad")), lwd = 2, lty = 1, col = "darkred")
# 2.3 ---- Testing differences among leaf types -------------------------------------------
# 2.3.1 - Parameter "A" ----
contrast(emmeans(M1, ~leaf.type, param = "A"), "pairwise")
contrast estimate SE df t.ratio p.value
Grass W - Grass E -30.9 8.84 226 -3.500 0.0031
Grass W - M. spinosum -223.2 15.92 226 -14.027 <.0001
Grass W - S. bracteolactus -240.3 16.86 226 -14.254 <.0001
Grass E - M. spinosum -192.3 13.45 226 -14.294 <.0001
Grass E - S. bracteolactus -209.4 14.22 226 -14.718 <.0001
M. spinosum - S. bracteolactus -17.0 11.71 226 -1.454 0.4671
P value adjustment: tukey method for comparing a family of 4 estimates
# 2.3.2 - Parameter "w" ----
contrast(emmeans(M1, ~leaf.type, param = "w"), "pairwise")
contrast estimate SE df t.ratio p.value
Grass W - Grass E 20.7 2.59 226 7.981 <.0001
Grass W - M. spinosum -31.7 7.75 226 -4.085 0.0004
Grass W - S. bracteolactus -27.0 8.07 226 -3.342 0.0053
Grass E - M. spinosum -52.4 6.62 226 -7.914 <.0001
Grass E - S. bracteolactus -47.7 6.84 226 -6.974 <.0001
M. spinosum - S. bracteolactus 4.7 6.72 226 0.699 0.8974
P value adjustment: tukey method for comparing a family of 4 estimates
Response function: lfmc = (A - w) / (1 + exp((m - time)/s))) + w
# 3.1.1 - Base nonlinear model ----
nl.fe.0 <- nls(lfmc ~ ((A - w) / (1 + exp((m - time) / s))) + w,
start = c(A = 15, m = 30, s = -17, w = 5),
data = lfmc
)
# here, time is the only predictor variable
nl.fe.0
Nonlinear regression model
model: lfmc ~ ((A - w)/(1 + exp((m - time)/s))) + w
data: lfmc
A m s w
194.97 29.66 -21.17 27.93
residual sum-of-squares: 1029759
Number of iterations to convergence: 8
Achieved convergence tolerance: 9.57e-07
b <- coef(nl.fe.0) # coefficients of the base model
b
A m s w
194.97252 29.66004 -21.17131 27.92798
# 3.1.2 - Leaf type fixed effect ----
nl.fe.1 <- nls(lfmc ~ ((A[leaf.type] - w[leaf.type]) / (1 + exp((m - time) / s))) + w[leaf.type],
start = list(A = rep(b[1], 4), m = 25, s = -10, w = rep(b[4], 4)),
data = lfmc
) # and fit is made using the base model's coefficients as start values
nl.fe.1
Nonlinear regression model
model: lfmc ~ ((A[leaf.type] - w[leaf.type])/(1 + exp((m - time)/s))) + w[leaf.type]
data: lfmc
A1 A2 A3 A4 m s w1 w2 w3 w4
56.726 92.483 298.835 318.207 30.431 -20.912 27.172 2.795 44.528 38.298
residual sum-of-squares: 67786
Number of iterations to convergence: 5
Achieved convergence tolerance: 5.958e-06
b1 <- coef(nl.fe.1) # coefficients of the model with leaf type and time as predictor variables
b1
A1 A2 A3 A4 m s w1 w2 w3 w4
56.725716 92.482865 298.834954 318.207067 30.430747 -20.911576 27.171667 2.795385 44.527958 38.297817
# 3.1.3 - Plot fixed effect ----
nl.fe.2 <- nls(lfmc ~ ((A[group] - w[group]) / (1 + exp((m - time) / s))) + w[group],
start = list(A = rep(c(b1[1], b1[2], b1[3], b1[4]), 3), m = 25, s = -10, w = rep(c(b1[7], b1[8], b1[9], b1[10]), 3)),
data = lfmc
) # the model is fit using "b1" as the start values
nl.fe.2
Nonlinear regression model
model: lfmc ~ ((A[group] - w[group])/(1 + exp((m - time)/s))) + w[group]
data: lfmc
A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 m s
53.683 54.282 61.906 111.562 314.131 333.831 77.465 298.161 321.433 87.553 279.437 292.443 30.413 -20.696
w1 w2 w3 w4 w5 w6 w7 w8 w9 w10 w11 w12
28.369 27.480 25.962 0.890 43.524 41.223 6.091 26.610 39.495 2.265 66.661 38.222
residual sum-of-squares: 53914
Number of iterations to convergence: 6
Achieved convergence tolerance: 7.249e-06
IC_tab(nl.fe.0, nl.fe.1, nl.fe.2)
# Residual checking:
par(mfrow = c(1, 1))
plot(nl.fe.2) # heterocedasticity in the residuals
hist(resid(nl.fe.2), main = "", xlab = "Residuals") # slightly skew distribution
qqnorm(resid(nl.fe.2))
qqline(resid(nl.fe.2))
This model assumes:
M2 <- nl.fe.2
summary(M2)
Formula: lfmc ~ ((A[group] - w[group])/(1 + exp((m - time)/s))) + w[group]
Parameters:
Estimate Std. Error t value Pr(>|t|)
A1 53.683 8.602 6.240 2.20e-09 ***
A2 54.282 8.636 6.285 1.72e-09 ***
A3 61.906 8.881 6.971 3.59e-11 ***
A4 111.561 13.291 8.394 5.65e-15 ***
A5 314.131 24.840 12.646 < 2e-16 ***
A6 333.831 26.640 12.531 < 2e-16 ***
A7 77.465 10.997 7.044 2.34e-11 ***
A8 298.161 24.917 11.966 < 2e-16 ***
A9 321.433 25.765 12.475 < 2e-16 ***
A10 87.553 11.744 7.455 2.01e-12 ***
A11 279.437 20.829 13.416 < 2e-16 ***
A12 292.443 24.178 12.095 < 2e-16 ***
m 30.413 2.896 10.504 < 2e-16 ***
s -20.695 3.606 -5.739 3.11e-08 ***
w1 28.369 6.537 4.340 2.17e-05 ***
w2 27.480 6.548 4.197 3.93e-05 ***
w3 25.962 6.633 3.914 0.000121 ***
w4 0.890 8.173 0.109 0.913389
w5 43.524 13.567 3.208 0.001535 **
w6 41.223 14.429 2.857 0.004686 **
w7 6.091 7.262 0.839 0.402467
w8 26.610 13.604 1.956 0.051725 .
w9 39.495 14.009 2.819 0.005252 **
w10 2.265 7.551 0.300 0.764455
w11 66.661 11.478 5.808 2.19e-08 ***
w12 38.222 13.031 2.933 0.003708 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 15.62 on 221 degrees of freedom
Number of iterations to convergence: 6
Achieved convergence tolerance: 7.249e-06
Response function: lfmc = β0 + β1xTime + β2xGE + β3xSM + β4xSS β5xGExTime + β6xSMxTime + β7xSSxTime
where GE, SM, and SS are dummy variables (0 or 1) created to indicate differences between the base level of “leaf type”, i.e., grasses from the W site [GW], and the remaining levels.
# 4.1.1 - Base linear mixed-effect model ----
l.me <- lme(lfmc ~ time * leaf.type, random = ~ 1 | plot, data = lfmc)
# Residual checking:
plot(l.me) # heterocedasticity in the residuals
plot(ACF(l.me, resType = "normalized", na.action = na.omit), alpha = 0.05, grid = TRUE) # temporal correlation is observed at lag 1
# 4.1.2 - Variance modelling ----
l.me.vm <- update(l.me, weights = varIdent(form = ~ 1 | leaf.type))
# Residual checking:
plot(l.me.vm) # variances there seem to be well modeled
AIC(l.me, l.me.vm) # modeling the variance improves the model fit
# 4.1.3 - Temporal correlation modelling ----
l.me.vm.tm <- update(l.me.vm, correlation = corARMA(p = 1, q = 1, form = ~1))
plot(ACF(l.me.vm.tm, resType = "normalized", na.action = na.omit), alpha = 0.05, grid = TRUE) # the temporal dependence is modeled
AIC(l.me.vm, l.me.vm.tm) # modeling the temporal correlation improves the model fit
plot(l.me.vm.tm)
plot(ACF(l.me.vm.tm, resType = "normalized", na.action = na.omit), alpha = 0.05, grid = TRUE)
hist(resid(l.me.vm.tm, type = "normalized"), main = "", xlab = "Residuals")
qqnorm(resid(l.me.vm.tm, type = "normalized"))
qqline(resid(l.me.vm.tm, type = "normalized"))
# 4.1.4 - Evaluation of the fixed-effects ----
l.me.ml <- update(l.me.vm.tm, method = "ML") # the model is fitted by Maximum likelihood
l.me.ml.1 <- update(l.me.ml, ~ . - time:leaf.type) # the model is reduced removing the interaction "time x leaf type"
IC_tab(l.me.ml, l.me.ml.1) # the AIC comparison suggests the interaction term to be important
M3 <- l.me.vm.tm
summary(M3)
Linear mixed-effects model fit by REML
Data: lfmc
AIC BIC logLik
1998.483 2050.63 -984.2416
Random effects:
Formula: ~1 | plot
(Intercept) Residual
StdDev: 3.750488 7.842661
Correlation Structure: ARMA(1,1)
Formula: ~1 | plot
Parameter estimate(s):
Phi1 Theta1
0.9076001 -0.7567060
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | leaf.type
Parameter estimates:
Grass W M. spinosum S. bracteolactus Grass E
1.000000 2.963571 3.123740 1.070897
Fixed effects: lfmc ~ time * leaf.type
Value Std.Error DF t-value p-value
(Intercept) 50.14877 3.478211 234 14.417977 0
time -0.27535 0.048114 234 -5.722819 0
leaf.typeGrass E 24.40614 4.908714 234 4.972002 0
leaf.typeM. spinosum 194.85292 8.539072 234 22.818980 0
leaf.typeS. bracteolactus 209.78398 8.838810 234 23.734415 0
time:leaf.typeGrass E -0.56367 0.071704 234 -7.861042 0
time:leaf.typeM. spinosum -2.02767 0.152372 234 -13.307359 0
time:leaf.typeS. bracteolactus -2.29338 0.160937 234 -14.250193 0
Correlation:
(Intr) time lf.tGE lf.M.s lf.S.b tm:.GE t:.M.s
time -0.557
leaf.typeGrass E -0.709 0.395
leaf.typeM. spinosum -0.407 0.227 0.654
leaf.typeS. bracteolactus -0.394 0.219 0.653 0.666
time:leaf.typeGrass E 0.374 -0.671 -0.592 -0.430 -0.429
time:leaf.typeM. spinosum 0.176 -0.316 -0.299 -0.742 -0.409 0.546
time:leaf.typeS. bracteolactus 0.166 -0.299 -0.318 -0.441 -0.743 0.562 0.567
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.3449313 -0.6438200 -0.1011963 0.6158869 3.3469739
Number of Observations: 247
Number of Groups: 6
M4 <- lm(lfmc ~ time * group, data = lfmc)
summary(M4)
Call:
lm(formula = lfmc ~ time * group, data = lfmc)
Residuals:
Min 1Q Median 3Q Max
-49.673 -7.042 0.466 7.002 66.332
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 48.79275 6.44309 7.573 9.60e-13 ***
time -0.24478 0.13305 -1.840 0.06714 .
groupGW plot 5 0.49817 9.11190 0.055 0.95645
groupGW plot 6 6.14984 9.11190 0.675 0.50042
groupGE plot 1 40.32675 9.49655 4.246 3.19e-05 ***
groupSM plot 1 212.68278 9.11190 23.341 < 2e-16 ***
groupSS plot 1 227.42468 9.11190 24.959 < 2e-16 ***
groupGE plot 2 14.39850 9.49655 1.516 0.13089
groupSM plot 2 195.54780 9.11190 21.461 < 2e-16 ***
groupSS plot 2 217.12857 9.11190 23.829 < 2e-16 ***
groupGE plot 3 21.68513 9.49655 2.283 0.02334 *
groupSM plot 3 189.79613 9.49655 19.986 < 2e-16 ***
groupSS plot 3 192.53490 9.49655 20.274 < 2e-16 ***
time:groupGW plot 5 -0.01905 0.18816 -0.101 0.91944
time:groupGW plot 6 -0.10231 0.18816 -0.544 0.58718
time:groupGE plot 1 -0.80129 0.19357 -4.139 4.94e-05 ***
time:groupSM plot 1 -2.36256 0.18816 -12.556 < 2e-16 ***
time:groupSS plot 1 -2.55767 0.18816 -13.593 < 2e-16 ***
time:groupGE plot 2 -0.43459 0.19357 -2.245 0.02574 *
time:groupSM plot 2 -2.34722 0.18816 -12.474 < 2e-16 ***
time:groupSS plot 2 -2.45552 0.18816 -13.050 < 2e-16 ***
time:groupGE plot 3 -0.56657 0.19357 -2.927 0.00378 **
time:groupSM plot 3 -1.82099 0.19357 -9.407 < 2e-16 ***
time:groupSS plot 3 -2.16847 0.19357 -11.202 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 16.51 on 223 degrees of freedom
Multiple R-squared: 0.9586, Adjusted R-squared: 0.9544
F-statistic: 224.7 on 23 and 223 DF, p-value: < 2.2e-16
# Residual checking:
par(mfrow = c(2, 2))
plot(M4) # a clear pattern in the residuals is observed (model assumptions are not met)
# 6 - NULL MODEL (M5) ######################################################################################
M5 <- lm(lfmc ~ 1, data = lfmc)
summary(M5)
Call:
lm(formula = lfmc ~ 1, data = lfmc)
Residuals:
Min 1Q Median 3Q Max
-88.75 -58.05 -31.45 52.66 231.06
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 95.645 4.919 19.45 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 77.3 on 246 degrees of freedom
IC_tab(M2, M3, M4, M5) # , weights = T, delta = TRUE, base=T, sort = TRUE) # the nonlinear mixed-effects model is clearly the best
IC_tab(M1)
par(mfrow = c(2, 2), cex = 0.65) # (Figure 5)
v1 <- c(1, 2, 3, 4)
v2 <- c("GW", "GE", "SM", "SS")
# Nonlinear mixed-effects model (M1)
rM1 <- resid(M1, type = "normalized")
fM1 <- fitted(M1)
plot(fM1, rM1, xlab = "Fitted LFMC (%)", ylab = "Residuals", ylim = c(-3, 3))
abline(a = 0, b = 0, lwd = 2)
boxplot(rM1 ~ lfmc$leaf.type, ylab = "Residuals", xlab = "Leaf type", xaxt = "n", ylim = c(-3, 3))
axis(side = 1, at = v1, labels = v2, las = 1)
plot(rM1 ~ lfmc$time, ylab = "Residuals", xlab = "Time (days)", ylim = c(-3, 3))
abline(a = 0, b = 0, lw = 2)
qqnorm(rM1, main = "", xlab = "Theoretical quantiles", ylab = "Sample quantiles")
qqline(rM1)
# Nonlinear fixed-effects model (M2)
par(mfrow = c(2, 2), cex = 0.65) # (Figure 5)
rM2 <- resid(M2)
fM2 <- fitted(M2)
plot(fM2, rM2, xlab = "Fitted LFMC (%)", ylab = "Residuals")
abline(a = 0, b = 0, lwd = 2)
boxplot(rM2 ~ lfmc$leaf.type, ylab = "Residuals", xlab = "Leaf type", xaxt = "n")
axis(side = 1, at = v1, labels = v2, las = 1)
plot(rM2 ~ lfmc$time, ylab = "Residuals", xlab = "Time (days)")
abline(a = 0, b = 0, lw = 2)
qqnorm(rM2, main = "", xlab = "Theoretical quantiles", ylab = "Sample quantiles")
qqline(rM2)
par(mfrow = c(2, 2), cex = 0.65)
# Linear mixed-effects model (M3)
rM3 <- resid(M3, type = "normalized")
fM3 <- fitted(M3)
plot(fM3, rM3, xlab = "Fitted LFMC (%)", ylab = "Residuals", ylim = c(-3, 3))
abline(a = 0, b = 0, lwd = 2)
boxplot(rM3 ~ lfmc$leaf.type, ylab = "Residuals", xlab = "Leaf type", xaxt = "n", ylim = c(-3, 3))
axis(side = 1, at = v1, labels = v2, las = 1)
plot(rM3 ~ lfmc$time, ylab = "Residuals", xlab = "Time (days)", ylim = c(-3, 3))
abline(a = 0, b = 0, lw = 2)
qqnorm(rM3, main = "", ylab = "Sample quantiles", xlab = "Theoretical quantiles")
qqline(rM3)
par(mfrow = c(2, 2), cex = 0.65)
# Clasical regression (M4)
rM4 <- resid(M4)
fM4 <- fitted(M4)
plot(fM4, rM4, xlab = "Fitted LFMC (%)", ylab = "Residuals")
abline(a = 0, b = 0, lwd = 2)
boxplot(rM4 ~ lfmc$leaf.type, ylab = "Residuals", xlab = "", xaxt = "n")
axis(side = 1, at = v1, labels = v2, las = 1)
plot(rM4 ~ lfmc$time, ylab = "Residuals", xlab = "Time (days)")
abline(a = 0, b = 0, lw = 2)
qqnorm(rM4, main = "", xlab = "Theoretical quantiles", ylab = "Sample quantiles")
qqline(rM4)